#### Distinguishing Exposure to Secondhand and Thirdhand Tobacco Smoke among U.S. Children Using Machine Learning: NHANES 2013−2016 ####
# https://pubs.acs.org/doi/10.1021/acs.est.2c08121?ref=PDF 原文链接
# 2023年3月30日08:45:11 进行复现论文paper数据
# 本文主要讨论区分二手烟和三手烟暴露区分开的问题  选取的数据是ehanes 2013-2014 2015-2016年数据
# 文章摘要说了关于
library(haven) #用于读取xpt
library(plyr) #用于数据处理包
library(dplyr) #用于数据处理包
library(arsenal) ##用于tableby
library(survey) #用于加权分析
library(stringr) #用于字符处理
library(gtsummary)
library(missRanger)
library(mltools)
library(caret)
library(randomForest)
library(sjPlot)
library(sjlabelled)
library(sjmisc)
library(ggplot2)
library(forestplot)


setwd("G:/BaiduNetdiskDownload")
#### 1 筛选数据   人口学数据 ####
# https://baike.baidu.com/item/%E5%8F%AF%E6%9B%BF%E5%AE%81/6702345
# 可替宁是尼古丁在人体内进行初级代谢后的主要产物——烟草中的尼古丁在体内经细胞色素氧化酶2A6(CYP2A6)代谢后的产物，
#主要存在于血液中，随着代谢过程从尿液排出。可替宁有促进神经系统兴奋作用，
#并在某些鼠类试验中反映出一定的抗炎、减轻肺水肿程度的作用。由于可替宁的半衰期较长(3~4d)且较稳定,因此成为测量吸烟者和被动吸烟者吸烟量的主要生物标志，
#一般情况下，多以血清中的可替宁浓度来评价。有研究成果显示，血浆中的可替宁浓度与血清中的可替宁浓度具有一致性，同样具有检测意义
# 方法 We fitted and tested random forest models, and the majority
#(76%) of children were classified in NEG, 16% were classified in TEG, and 8% were classified in
# MEG  文章使用了随机检测森林模型 76%的人是在meg 16% 的在TEG 最后8% 在meg
# 参数解读
# NEG	no/minimal tobacco smoke exposure	极少接触暴露在烟草
# TEG	predominant THS exposure	主要三手烟暴露
# MEG	mixed SHS and THS exposure	混杂的二手烟和三手烟暴露
# 2023年4月28日15:25:42 
# 参与者是4485名不吸烟的3-17岁青少年，来自国家健康与营养部2013-2016年考试调查
# 我们拟合并测试了随机森林模型（76%）的儿童被分类为NEG，16%被分类为TEG，8%被分类为MEG
# 结论 最终基于报告 生物标志物最终分类模型预测准确率95% 模型在NEG准确率 100% TEG 准确率88% MEG准确率 71%
#最重要预测因素有 家庭吸烟者数量(number of household smokers) ,血清可替宁 (serum cotinine) ,血清羟基可替宁(serum hydroxycotinine)
# , 尿4(urinary 4-(methylnitrosamino)-1-(3-pyridyl)-1-butanol (NNAL))

# 在缺乏经过验证的特定于 THS 的生物标志物的情况下，全面的生物标志物和烟草烟雾暴露问卷数据可以高精度地区分暴露于 SHS 和 ths 的儿童。。

# 2023年4月28日16:02:09 
#总尼古丁当量（TNEs）是 被认为是每天估计的黄金标准生物标志物
# 2023年5月3日15:45:56 
# 刚把全文翻译了进行 目前还不知道y变量是什么 (ಥ_ಥ) 
# 在摘要最后一句说了综合的生物标志物和烟草暴露的问卷调查数据可以提高儿童区分二手烟和三手烟精确度

##### 1.1  通过表格1 可以看到涉及 年龄 种族 性别 学历 等信息 #####

demo_h <- read_xpt("NHANES/2013-2014/Demographics/demo_h.xpt")
demo_i <- read_xpt("NHANES/2015-2016/Demographics/demo_i.xpt")
#合并两个cycle 数据
demo_data_all <- dplyr::bind_rows(list(demo_h,demo_i))
demo_data <- demo_data_all[,c('SEQN', 'RIDAGEYR', 'RIAGENDR','RIDRETH1', 'RIDRETH3', 'DMDEDUC2','DMDHHSIZ', 'INDFMPIR' ,'DMDHHSZA','DMDHHSZB' )]

# 进行测试数据
# demo_datalower12 <-subset.data.frame(demo_data,RIDAGEYR<=11)
# demo_dataover12 <-subset.data.frame(demo_data,RIDAGEYR>11)
# table(demo_datalower12$SIAPROXY)
# 
# table(demo_dataover12$SIAPROXY)


# We delimited our analysis to 4485 nonsmoking children who
# were aged 3−17 years and had serum cotinine results available for analysis

# 文章中说最后分析的结果限定为4485个没有吸烟的3-17岁儿童且有血清可替宁检测到的数据
# 排除过程 
# NHANES did not measure the biomarkers included in this study on
#participants <3 years old, and therefore, the younger group
#was excluded from analyses (n = 2758).  Children missing data
#on household smokers and home TSE were also excluded prior
#to analyses (n = 70).  Additionally, children ≥12 years old (n =3038) were asked if they smoked any tobacco product (i.e.,cigarettes, electronic cigarettes 
#[e-cigarettes], cigars, littlecigars/cigarillos, hookah, waterpipes) in the past 5 days or ifthey smoked cigarettes in the past 30-days and were excluded if
#they reported past 5-day tobacco use or past 30-day cigarette
#smoking (n = 94).

##### 1.2 SMQ-吸烟数据提取 #####
### 1.1 提取 Component文件 
# 这里需要提取 Cigarette Use ,Smoking - Household Smokers,Smoking - Secondhand Smoke Exposure
# Questionnaire	hoq_h

# 看到有个监护人,但是我没找到怎么搞 就是全部数据

smq_h <- read_xpt("NHANES/2013-2014/Questionnaire/smq_h.xpt")
smq_i <- read_xpt("NHANES/2015-2016/Questionnaire/smq_i.xpt")
smqfam_h <- read_xpt("NHANES/2013-2014/Questionnaire/smqfam_h.xpt")
smqfam_i <- read_xpt("NHANES/2015-2016/Questionnaire/smqfam_i.xpt")
smqshs_h <- read_xpt("NHANES/2013-2014/Questionnaire/smqshs_h.xpt")
smqshs_i <- read_xpt("NHANES/2015-2016/Questionnaire/smqshs_i.xpt")
# 询问家里房子多少 
hoq_h <- read_xpt("NHANES/2013-2014/Questionnaire/hoq_h.xpt")
hoq_i <- read_xpt("NHANES/2015-2016/Questionnaire/hoq_i.xpt")

# SMQRTU
smqrtu_h <- read_xpt("NHANES/2013-2014/Questionnaire/smqrtu_h.xpt")
smqrtu_i <- read_xpt("NHANES/2015-2016/Questionnaire/smqrtu_i.xpt")

smq_all <- dplyr::bind_rows(list(smq_h,smq_i))
smqfam_all <- dplyr::bind_rows(list(smqfam_h,smqfam_i))
smqshs_all <- dplyr::bind_rows(list(smqshs_h,smqshs_i))
hoq_all <- dplyr::bind_rows(list(hoq_h,hoq_i))
smqrtu_all <- dplyr::bind_rows(list(smqrtu_h,smqrtu_i))


smq_data_all <- plyr::join_all(list(smq_all,smqfam_all,smqshs_all,hoq_all,smqrtu_all),
                         by='SEQN',type='full')

#table(smq_data_all$SMAQUEX)

#outputlower11 <-subset.data.frame(smq_data_all)

# 进行验证数量
#SMAQUEX 是问卷模式

smq_data <- smq_data_all[,c('SEQN', 'SMD460', 'SMD470', 'SMD641', 'SMD480','HOD050','HOQ065','SMDANY','SMQ860','SMQ862','SMQ870','SMQ872','SMQ874','SMQ876','SMQ878','SMQ880')]

#### 1.3 lbdcotlc 可替宁提取####
### 1.1 提取 Component文件 Y变量
#
cot_h <- read_xpt("NHANES/2013-2014/Laboratory/cot_h.xpt")
cot_i <- read_xpt("NHANES/2015-2016/Laboratory/cot_i.xpt")
# 提起实验室部分数据 ->
pah_h <- read_xpt("NHANES/2013-2014/Laboratory/pah_h.xpt")
pah_i <- read_xpt("NHANES/2015-2016/Laboratory/pah_i.xpt")

ucot_h <- read_xpt("NHANES/2013-2014/Laboratory/ucot_h.xpt")
ucot_i <- read_xpt("NHANES/2015-2016/Laboratory/ucot_i.xpt")
tsna_h <- read_xpt("NHANES/2013-2014/Laboratory/tsna_h.xpt")

# UVOC
uvoc_h <- read_xpt("NHANES/2013-2014/Laboratory/uvoc_h.xpt")
uvoc_i <- read_xpt("NHANES/2015-2016/Laboratory/uvoc_i.xpt")



cot_data_all <- dplyr::bind_rows(list(cot_h,cot_i))
pah_data_all <- dplyr::bind_rows(list(pah_h,pah_i))
ucot_data_all <- dplyr::bind_rows(list(ucot_h,ucot_i))
uvoc_data_all <- dplyr::bind_rows(list(uvoc_h,uvoc_i))
#colnames(cot_data_all)
# 2023年5月20日14:43:21 进行提取需要的数据

Laboratory_data_all <- plyr::join_all(list(cot_data_all,pah_data_all,ucot_data_all,tsna_h,uvoc_data_all),
                               by='SEQN',type='full')

# 血清可替宁->lbxcot
# 羟可替宁，血清 (ng/mL)->  lbxhct
# 总羟基可替宁，尿液 (ng/mL) ->  urxhctt
# 总可替宁，尿液 (ng/mL) ->  urxcott
# NNAL，尿液 (ng/mL) -> urxnal --tsna_h.xpt
# 2-羟基芴 (ng/L) -> urxp04
# 3-羟基芴 (ng/L) -> urxp03
# N-乙酰基-S-(2-氰乙基)-L-半胱氨酸(ng/mL) -> urxcym -- 
cot_data <- Laboratory_data_all[,c('SEQN', 'LBXCOT','LBXHCT','URXHCTT','URXCOTT','URXNAL','URXP04','URXP03','URXCYM' )]

#### 组合全部数据 先这样 后续补充如果缺少的
output <- plyr::join_all(list(demo_data, smq_data_all, cot_data),
                         by='SEQN',type='full')

outputdata <- output[,c('SEQN','RIDAGEYR','DMDEDUC2','INDFMPIR', 'SMD641', 'DMDHHSIZ', 'SMDANY', 'DMDHHSZA','DMDHHSZB','HOQ065','LBXCOT','LBXHCT', 'URXHCTT','URXCOTT','URXNAL','URXP04','URXP03','URXCYM','RIDRETH1', 'HOD050',  'RIAGENDR','SMAQUEX','SMD460','SMD470', 'SMD480','HOD050','HOQ065','SMQ860','SMQ862','SMQ870','SMQ872','SMQ874','SMQ876','SMQ878','SMQ880')]

# write.csv(outputdata,'D:/Users/GCR/Desktop/需要复现文章三篇/outputdata.csv')
#outputdata<-na.omit(outputdata)
#table(outputdata$RIDAGEYR)
#table(outputdata$SMAQUEX)
# 2023年5月11日14:41:32 搞出来了 就是这个方法 上面的数据就是代理者访问 里面的数据都是3 也就是根据nhanes的解释 3 是代理0-11岁回答的 里面数据分布
#   0   1   2   3   4   5   6   7   8   9  10  11 
#  405 262 278 220 230 208 262 233 217 213 223 223 
#dim(outputdata) #0-11 岁 儿童数据
#output18SMAQUEX322 <-subset.data.frame(outputdata,SMAQUEX==3&RIDAGEYR>=3&RIDAGEYR<=11)
#output18SMAQUEX3222 <-subset.data.frame(outputdata,SMAQUEX==2)
#a3 <- output18SMAQUEX322[,'SEQN']
#a2 <- output18SMAQUEX3222[,'SEQN']
# 这两个没有交集 也就是不是通过 SMAQUEX 进行区分的人 我以为同一个seqn 通过不同类型的SMAQUEX 区分 但是不是
#aa <- intersect(a3,a2)
# 下一步处理12-17 岁儿童
# outputdataage1217 <-  subset.data.frame(outputdata,RIDAGEYR>=12&RIDAGEYR<=17)

# table(outputdataage1217$SMAQUEX)
# table(outputdataage1217$RIDAGEYR)
# size1217 <- dim(outputdataage1217)
# size311 <- dim(output18SMAQUEX322)
# 开始测试有完整数据的 2023年5月11日15:26:15

#child age, M (SE) 
# 直接过滤3-17岁就完了
# child sex
outputdata$Sex <- ifelse(outputdata$RIAGENDR==1,'male','female')
# child race/ethnicity
race <- recode_factor(outputdata$RIDRETH1, 
                      `1` = 'Mexican American',
                      `2` = 'Hispanic',
                      `3` = 'Non-Hispanic White',
                      `4` = 'Non-Hispanic Black',
                      `5` = 'Other Race - Including Multi-Racial'
                 
)
outputdata$race <- race
# caregiver age, M (SE)
# 直接采用大于17岁完事

# caregiver sex 
# 与上面类似

#caregiver education level
# table(outputdata$DMDEDUC2)
edu <- recode_factor(outputdata$DMDEDUC2, 
                      `1` = '≤high school graduate/equivalent',
                      `2` = '≤high school graduate/equivalent',
                      `3` = '≤high school graduate/equivalent',
                      `4` = 'some college',
                      `5` = '≥college graduate'
                      
)
outputdata$caregivereducationlevel <- edu
# FPL INDFMPIR
# aaa<-which(is.na(outputdata$INDFMPIR))
# aaa
# outputdata[7922,'FPL']
outputdata$FPL <- ifelse(outputdata$INDFMPIR<=1.85,'<185%',ifelse((outputdata$INDFMPIR>=1.85&outputdata$INDFMPIR<=3.49) , "185−349%" ,'≥350%'))
# 添加NA的是未指定
outputdata$FPL[which(is.na(outputdata$FPL))] <- 'unspecified'

# household members, 家庭成员数量均数

# DMDHHSIZ

# child household members 家庭儿童数量
outputdata$allchild <- outputdata$DMDHHSZA+outputdata$DMDHHSZB

chindData <- recode_factor(outputdata$allchild, 
                     `1` = '1 child',
                     `2` = '2 children',
                     `3` = '≥3 children',
                     `4` = '≥3 children',
                     `5` = '≥3 children',
                     `6` = '≥3 children',
                     `7` = '≥3 children',
                     `8` = '≥3 children',
                     `9` = '≥3 children',
                     `10` = '≥3 children',
                 
                     
)
outputdata$childhouseholdmembers <- chindData
# family home ownership status
# Questionnaire	HOQ	hod050  HOD050 - Number of rooms in home  HOQ065 - Home owned, bought, rented, other

familyhomeownershipstatus <- recode_factor(outputdata$HOQ065, 
                           `1` = 'owns home/being bought',
                           `2` = 'rents home/other arrangement',
                           `3` = 'rents home/other arrangement'
                           
)
outputdata$familyhomeownershipstatus <- familyhomeownershipstatus
# str(outputdata)


# 原文提示在屋子里面吸烟的 对应的原文就是SMD460 这个参数 现在把 SMD460  SMAQUEX SEQN  RIDAGEYR 拿出来看下
# 通过检查年龄 问卷模式 室内吸烟的情况,初步判断出 监护人暂且不管了,找不到,进行数据后续操作
# 首先根据文章问卷调查把儿童分为三组 
# https://wwwn.cdc.gov/Nchs/Nhanes/2013-2014/SMQFAM_H.htm
#2013-2014/SMQFAM_H.htm  首先第一步 Smoking - Household Smokers 模块 初步筛选出人群 解决这块 Household Smokers and Home Tse 段落问题
# 进行分组 NEG TEG MEG 
# 首先是家里面没有吸烟的 初步分为NEG 
# SMD480 过去一周几天有人在这里吸烟

# 下面这个是错误的 需要多级测试没问题才可以使用
#PEPGROUP <- recode_factor(outputdata$SMD460, 
#                                           `0` = 'NEG'
#)
#outputdata$PEPGROUP <- PEPGROUP

#SHS Exposure in the Past 7 Days. To further distinguish
#children into the NEG, TEG, and MEG based on exposure in
#locations other than the home, we evaluated SHS exposure in
#various environments in the past 7 days via four items from the
#questionnaire’s “Smoking − SHS Exposure” section. For each
#SHS exposure item, a preliminary filter yes/no question was
#asked to determine whether the child spent time in the
#particular environment during the past 7 days (i.e., restaurant,car, another home, and any other indoor area). If the child was
#in an environment (e.g., rode in a car), then a follow-up yes/no
#question was asked about whether someone else smoked
#tobacco products while the participant was in that respective
#environment in the past 7 days. Participants with a “yes”
#response were considered exposed to SHS in the environment

# https://wwwn.cdc.gov/Nchs/Nhanes/2013-2014/SMQSHS_H.htm
# 2013-2014/SMQSHS_H.htm 这个模块可以 解决 SHS Exposure in the Past 7 Days 这块内容的需要的数据

# 暴露于SHS SHS Exposure in the Past 7 Days. #table2 部分采用数据来源
SMQ860 <- ifelse(outputdata$SMQ860==1&outputdata$SMQ862==1,'y','n')
SMQ870 <- ifelse(outputdata$SMQ870==1&outputdata$SMQ872==1,'y','n')
SMQ874 <- ifelse(outputdata$SMQ874==1&outputdata$SMQ876==1,'y','n')
SMQ878 <- ifelse(outputdata$SMQ878==1&outputdata$SMQ880==1,'y','n')

outputdata$cartse <- as.factor(SMQ870)
outputdata$otherhomeTSE <- as.factor(SMQ874)
outputdata$restaurantTSE <- as.factor(SMQ860)
outputdata$otherindoorareaTSE <- as.factor(SMQ878)

# 设置表2 SMD460 smoker判断
table(outputdata$SMD460)
smokersize <- recode_factor(outputdata$SMD460, 
                           `0` = '0 smokers',
                           `1` = '1 smokers',
                           `2` = '≥2 smokers',
                           `3` = '≥2 smokers'
                           
                        
)
outputdata$smokersize <- smokersize
SHSenvironment <- ifelse(SMQ860=='y'| SMQ870=='y'| SMQ874=='y'| SMQ878=='y',1,0)
outputdata$SHSenvironment <-SHSenvironment
# 如果生活在不与吸烟居住一起 且最近一周shs 没有的 则NEG 有的是MEG 
# 如果与吸烟住一起 然后判断是不是在家里面吸烟 no的 进行判断二手烟暴露 no的 TEG 是的 是MEG  如果 家里面吸烟判断二手烟暴露 no的是meg 是的是meg
# lives with smoker ->   SMD460  smoked inside home -> SMD470  

outputdata$smd460Flag <- ifelse(outputdata$SMD460>=1&outputdata$SMD460<=3,1,ifelse(outputdata$SMD460==0,0,''))
outputdata$smd470Flag <- ifelse(outputdata$SMD470>=1&outputdata$SMD470<=3,1,ifelse(outputdata$SMD470==0,0,''))
# 图1 流程图转换如下 2023年5月15日19:15:27 
#smd460Flag #SHSenvironment  smd470Flag   #SHSenvironment  result
# 0               0                                         NEG
# 0               1                                         MEG 
# 1                            0                 0          TEG
# 1                            0                 1          MEG
# 1                            1                 0          MEG
# 1                            1                 1          MEG

outputdata$pgroup1<- ifelse(outputdata$smd460Flag==0&outputdata$SHSenvironment==0,'NEG',"")
outputdata$pgroup2<- ifelse(outputdata$smd460Flag==0&outputdata$SHSenvironment==1,'MEG',"")
outputdata$pgroup3<- ifelse(outputdata$smd460Flag==1&outputdata$smd470Flag==0&outputdata$SHSenvironment==0,'TEG',"")
outputdata$pgroup4<- ifelse(outputdata$smd460Flag==1&outputdata$smd470Flag==0&outputdata$SHSenvironment==1,'MEG',"")
outputdata$pgroup5<- ifelse(outputdata$smd460Flag==1&outputdata$smd470Flag==1&outputdata$SHSenvironment==0,'MEG',"")
outputdata$pgroup6<- ifelse(outputdata$smd460Flag==1&outputdata$smd470Flag==1&outputdata$SHSenvironment==1,'MEG',"")
outputdata$allGroup<- paste(outputdata$pgroup1,outputdata$pgroup2,outputdata$pgroup3,outputdata$pgroup4,outputdata$pgroup5,outputdata$pgroup6,collapse = NULL)
outputdata$allGroup <- str_replace_all(outputdata$allGroup,' ','')
outputdata$allGroup <- str_replace_all(outputdata$allGroup,'NA','')

# 先复现 儿童的 在复现监护人
agerange317<-subset.data.frame(outputdata,RIDAGEYR>=3&RIDAGEYR<=17)
dim(agerange317) #6116
# 删除不含有血清可替宁的人

agerange317dropna <- subset.data.frame(agerange317,!is.na(LBXCOT))
dim(agerange317dropna) # LBXCOT 4649
# 删除最近五天吸烟的人 或者过去三十天吸烟的人
# 2013-2014/SMQRTU_H.htm
# SMDANY - Used any tobacco product last 5 days?  1 删除
# SMD641 # days smoked cigs during past 30 days -> On how many of the past 30 days did {you/SP} smoke a cigarette? 1-30 删除
# table(agerange317dropna$SMD641)
# table(agerange317dropna$SMDANY)
indexsmd641del <- which(agerange317dropna$SMD641>=1)
# length(indexsmd641del)
agerange317dropna <- agerange317dropna[-indexsmd641del,]
# dim(agerange317dropna)
indexSMDANYdel <- which((agerange317dropna$SMDANY==1))
# length(indexSMDANYdel)
agerange317dropna <- agerange317dropna[-indexSMDANYdel,]
# dim(agerange317dropna) 

#删除没有TSE 我认为是就是 allGroup 是空额度
table(agerange317dropna$allGroup)

indexgroupdel <- which(agerange317dropna$allGroup=='')
# length(indexgroupdel)
agerange317dropna <- agerange317dropna[-indexgroupdel,]

write.csv(agerange317dropna,'D:/Users/GCR/Desktop/需要复现文章三篇/outputdata.csv')
dim(agerange317dropna) #4302    2023年5月16日14:10:31 还是有点差距 先不管了 实在是找不出了
colnames(agerange317dropna)
# agerange317dropna0.03 <- subset.data.frame(agerange317,LBXCOT>=0.03)
# table(agerange317dropna$LBXHCT)
# 测试表格1  暂时不知道怎么加95% 
str(agerange317dropna)
table(agerange317dropna$allGroup)
agerange317dropna$allGroup <-as.factor(agerange317dropna$allGroup)

# tblchild1<- tbl_summary(agerange317dropna,
#             missing = 'n',
#             # 通过属性分组
#             by = allGroup,
#             # 进行设置显示维度
#             include = c(RIDAGEYR,Sex,race,FPL,DMDHHSIZ,childhouseholdmembers,familyhomeownershipstatus,HOD050),
#             # 批量给变量起别名
#             label = list(RIDAGEYR ~ "child age, M (SE)",
#                          Sex ~ "child sex",
#                          race ~ "child race/ethnicity",
#                          FPL ~ "FPL",
#                          DMDHHSIZ ~ "no. household members, M (SE)",
#                          childhouseholdmembers ~ "no. child household members",
#                          familyhomeownershipstatus ~ "family home ownership status",
#                          HOD050 ~ "no. household rooms, M (SE)"
#                          ),
# 
#             statistic = list(
#               # 分别对应数值型和分类变量 后是需要显示表达式
#               all_continuous() ~ "{mean}",
#               all_categorical() ~ "{n} ({p} ) "
#             ),
#             digits = all_continuous() ~ 3
#             )%>%
#   #add_n() %>% # 添加非NA观测值个数
#   # 这里p值不能用p 根据原文说的是加权得来的 需要用 tbl_svysummay进行 这里除了p值 其他是正确的
#   # add_p() %>% # 添加P值
#   add_overall() %>%
#   #add_significance_stars(hide_se = FALSE,hide_p = TRUE)%>%
#   #modify_header(std.error = "**SE**")
#   add_ci() %>%
#   modify_column_merge(pattern = "{stat_0} ({ci_stat_0})")%>%
#   modify_column_merge(pattern = "{stat_1} ({ci_stat_1})")%>%
#   modify_column_merge(pattern = "{stat_2} ({ci_stat_2})")%>%
#   modify_column_merge(pattern = "{stat_3} ({ci_stat_3})")%>%
#   #modify_spanning_header(c("stat_2", "stat_3","stat_4") ~ "**TSE group**")%>%
#   modify_header(label = "**characteristic**",all_stat_cols() ~ "**{level}**<br>N = {n} <br> n(% [95% CI]")
# 
# tblchild1
  

# show_header_names(tblchild1)
# 表格2 和上面一样 暂时不知道怎么加95%
# 进行复现表格2 部分数据进行转换
# number of household smokers -> smokersize 




agerange317dropna$cartse <-as.factor(agerange317dropna$cartse)
agerange317dropna$otherhomeTSE <-as.factor(agerange317dropna$otherhomeTSE)
agerange317dropna$restaurantTSE <-as.factor(agerange317dropna$restaurantTSE)
agerange317dropna$otherindoorareaTSE <-as.factor(agerange317dropna$otherindoorareaTSE)

table2data<- agerange317dropna[,c('smokersize','cartse','otherhomeTSE','restaurantTSE','otherindoorareaTSE','allGroup')]

str(table2data)

tblchild2<-tbl_summary(table2data,
                       # 通过属性分组
                       by = allGroup,
                       # 进行设置显示维度
                       include = c(smokersize,cartse,otherhomeTSE,restaurantTSE,otherindoorareaTSE),
                       # 批量给变量起别名
                       label = list(smokersize ~ "number of household smokers",
                                    cartse ~ "car TSE",
                                    otherhomeTSE ~ "other home TSE",
                                    restaurantTSE ~ "restaurant TSE",
                                    otherindoorareaTSE ~ "other indoor area TSE"
                                    
                       ),
                      
                       
                       statistic = list(
                         # 分别对应数值型和分类变量 后是需要显示表达式
                         all_continuous() ~ "{mean}±{sd}",
                         all_categorical() ~ "{n} ({p}%)"
                       ),
                       digits = all_continuous() ~ 3
)%>%
  #add_n() %>% # 添加非NA观测值个数
  # 这里p值不能用p 根据原文说的是加权得来的 需要用 tbl_svysummay进行 这里除了p值 其他是正确的
  # add_p() %>% # 添加P值
  add_overall() %>%
  add_ci(pattern = "{stat} ({ci})") %>%
  #modify_spanning_header(c("stat_2", "stat_3","stat_4") ~ "**TSE group**")%>%
  modify_header(label = "**TSE variable**")

tblchild2
# 表格3  类似 先不管了  不行 后面随机森林用到这块参数
# 表格三主要是生物测量数据
# 血清尼古丁代谢物 
#1  血清可替宁 serum cotinine ->  lbxcot
#2 血清羟可替宁  serum hydroxycotinine -> lbxhct
#3 尿tne2 >-  urinary TNE2 ->  urxcott/176.2151 +  urxhctt/192.2145 
#4  NNAL (pg/mL) NNAL，尿液 (ng/mL) 
#5 NNAL/TNE2 -> #4/ #3
#6 2-羟基芴 (ng/L) -> 2-hydroxyfluorene (ng/L) -> urxp04
#7 3-羟基芴 (ng/L) -> 3-hydroxyfluorene (ng/L) -> urxp03
#8 2-hydroxyfluorene/TNE2 -> #6 / #3 
#9 3-hydroxyfluorene/TNE2 -> #7 / #3 
# 10 N-乙酰基-S-(2-氰乙基)-L-半胱氨酸(ng/mL) -> urxcym 2CyEMA
# 11 N-乙酰基-S-(2-氰乙基)-L-半胱氨酸(ng/mL) /TNE2 -> #10 /  #3 
 # 拿过来的数据
# 血清可替宁->lbxcot
# 羟可替宁，血清 (ng/mL)->  lbxhct
# 总羟基可替宁，尿液 (ng/mL) ->  urxhctt
# 总可替宁，尿液 (ng/mL) ->  urxcott
# NNAL，尿液 (ng/mL) -> urxnal --tsna_h.xpt
# 2-羟基芴 (ng/L) -> urxp04
# 3-羟基芴 (ng/L) -> urxp03
# N-乙酰基-S-(2-氰乙基)-L-半胱氨酸(ng/mL) -> urxcym --  2CyEMA
colnames(agerange317dropna)
#urinary TNE2
agerange317dropna$urinaryTNE2 <- (agerange317dropna$URXCOTT /176.2151) +(agerange317dropna$URXHCTT /192.2145  )
#NNAL/TNE2 
agerange317dropna$NNAL_TNE2 <- (agerange317dropna$URXNAL /agerange317dropna$urinaryTNE2)
#2-hydroxyfluorene (ng/L) 
agerange317dropna$twohydroxyfluorene<- agerange317dropna$URXP04 
#3-hydroxyfluorene (ng/L)
agerange317dropna$threehydroxyfluorene<- agerange317dropna$URXP03 
#2-hydroxyfluorene/TNE2 
agerange317dropna$twohydroxyfluorene_tne2<- agerange317dropna$twohydroxyfluorene /agerange317dropna$urinaryTNE2
#3-hydroxyfluorene/TNE2 
agerange317dropna$threehydroxyfluorene_tne2<- agerange317dropna$threehydroxyfluorene /agerange317dropna$urinaryTNE2
# 2CyEMA (ng/mL) 
agerange317dropna$twoCyEMA<- agerange317dropna$URXCYM 
# 2CyEMA/TNE2
agerange317dropna$twoCyEMA_TNE2<- agerange317dropna$twoCyEMA /agerange317dropna$urinaryTNE2
colnames(agerange317dropna)
# 血清可替宁->lbxcot
# 羟可替宁，血清 (ng/mL)->  lbxhct
# 总羟基可替宁，尿液 (ng/mL) ->  urxhctt
# 总可替宁，尿液 (ng/mL) ->  urxcott
# NNAL，尿液 (ng/mL) -> urxnal --tsna_h.xpt
# 2-羟基芴 (ng/L) -> urxp04
# 3-羟基芴 (ng/L) -> urxp03
# N-乙酰基-S-(2-氰乙基)-L-半胱氨酸(ng/mL) -> urxcym --  2CyEMA
# tblchild3<- tbl_summary(agerange317dropna,
#             missing = 'n',
#             # 通过属性分组
#             by = allGroup,
#             # 进行设置显示维度
#             include = c(LBXCOT,LBXHCT,urinaryTNE2,URXNAL,NNAL_TNE2,twohydroxyfluorene,threehydroxyfluorene,twohydroxyfluorene_tne2,threehydroxyfluorene_tne2,twoCyEMA,twoCyEMA_TNE2),
#             # 批量给变量起别名
#             label = list(
#               LBXCOT ~ "serum cotinine (ng/mL)",
#               LBXHCT ~ "serum hydroxycotinine (ng/mL)",
#               urinaryTNE2 ~ "urinary TNE2 (nmol/mL)",
#               URXNAL ~ "NNAL (pg/mL)",
#               NNAL_TNE2 ~ "NNAL/TNE2",
#               twohydroxyfluorene ~ "2-hydroxyfluorene (ng/L)",
#               threehydroxyfluorene ~ "3-hydroxyfluorene (ng/L)",
#               twohydroxyfluorene_tne2 ~ "2-hydroxyfluorene/TNE2",
#               threehydroxyfluorene_tne2 ~ "3-hydroxyfluorene/TNE2",
#               twoCyEMA ~ "2CyEMA (ng/mL)",
#               twoCyEMA_TNE2 ~ "2CyEMA/TNE2"),
# 
#             statistic = list(
#               # 分别对应数值型和分类变量 后是需要显示表达式
#               all_continuous() ~ "{mean}"
#              # all_categorical() ~ "{n} ({p}) "
#             ),
#             digits = all_continuous() ~ 2
#             )%>%
#   add_n(statistic = "{n}({p})") %>% # 添加非NA观测值个数
#   # 这里p值不能用p 根据原文说的是加权得来的 需要用 tbl_svysummay进行 这里除了p值 其他是正确的
#   # add_p() %>% # 添加P值
#   add_overall() %>%
#   #add_significance_stars(hide_se = FALSE,hide_p = TRUE)%>%
#   #modify_header(std.error = "**SE**")
#   add_ci() %>%
#   # modify_column_merge(pattern = "{stat_0} ({ci_stat_0})")%>%
#   modify_column_merge(pattern = "{stat_1} ({ci_stat_1})")%>%
#   modify_column_merge(pattern = "{stat_2} ({ci_stat_2})")%>%
#   modify_column_merge(pattern = "{stat_3} ({ci_stat_3})")%>%
#   modify_spanning_header(c("stat_1", "stat_2","stat_3") ~ "**TSE group membership**")%>%
#   modify_header(label = "**biomarker variable**",all_stat_cols() ~ "**{level}**<br>GeoM (95% CI) ")%>% add_overall()
# 
# 
# show_header_names(tblchild3)
# tblchild3
#agerange317dropna <- missRanger(agerange317dropna, pmm.k = 3, num.trees = 100)
#colnames(agerange317dropna)

# 主要表格4  随机森林模型 
# 还需要学习gini impurity 基尼系数 
# 16个重要变量进行测试
#  number of house smoker -> smokersize
#  serum cotinine ->LBXCOT
#  serum Hydroxycotinine ->LBXHCT
#  urinary nnal -> URXNAL
#  3-hydroxyfluorene/TNE2  -> threehydroxyfluorene_tne2
#  urinary TNE2 -> urinaryTNE2
#  2-hydroxyfluorene/TNE2  -> twohydroxyfluorene_tne2
#  2CyEMA/TNE2 -> twoCyEMA_TNE2
#  NNAL/TNE2  ->NNAL_TNE2
#  2-hydroxyfluorene -> twohydroxyfluorene
#  car tse ->cartse
#  3-hydroxyfluorene -> threehydroxyfluorene
#  2CyEMA  -> twoCyEMA
#  other home tse -> otherhomeTSE
#  otherindoorareaTSE -> otherindoorareaTSE
#  restaurantTSE -> restaurantTSE

colnames(agerange317dropna)
finalModel <- agerange317dropna[,c( 'smokersize','LBXCOT','LBXHCT','URXNAL','threehydroxyfluorene_tne2','urinaryTNE2',
                                   'twohydroxyfluorene_tne2','twoCyEMA_TNE2','NNAL_TNE2','twohydroxyfluorene','cartse','threehydroxyfluorene',
                                   'twoCyEMA','otherhomeTSE','otherindoorareaTSE','restaurantTSE','allGroup')]

# 文章说的提取
finalModel <- missRanger(finalModel, pmm.k = 2, num.trees = 500)
# table(finalModel$allGroup)
# dim(finalModel)
# 根据数据进行分析 拆分数据 80 训练 20 验证
# set.seed(40)
# trains <- createDataPartition(y = finalModel$allGroup,p = 0.8,list = FALSE)
# trainsdata <- finalModel[trains,]
# testdata <- finalModel[-trains,]
# 这块说的是下面六个模型进行拿着20%的数据进行验证 拿到结果进型六个额外模型验证 直接验证吧 没时间了
# all reported TSE and biomarker variables -> finalModel
# all reported TSE variables -> finalModel1
# reported TSE variables excluding number of household smokers -> finalMode2
# 所有的生物标志物排除了比值的
# all biomarker variables excluding ratios -> finalModel3
# biomarker and biomarker ratio variables excluding serum cotinine, serum hydroxycotinine, urinary TNE2, and urinary NNAL -> finalModel4
# all biomarker and biomarker ratio variables -> finalModel5
# all biomarker ratio variables -> finalModel6
# 首先生成这六额外个模型

finalModel1 <- finalModel[,c('otherhomeTSE','otherindoorareaTSE','restaurantTSE','cartse','smokersize','allGroup')]
finalModel2 <- finalModel[,c('otherhomeTSE','otherindoorareaTSE','restaurantTSE','cartse','allGroup')]
finalModel3 <- finalModel[,c('LBXCOT','LBXHCT','urinaryTNE2','URXNAL','twohydroxyfluorene','threehydroxyfluorene', 'twoCyEMA', 'allGroup')]
finalModel4 <- finalModel[,c('twohydroxyfluorene','threehydroxyfluorene', 'twoCyEMA','NNAL_TNE2','threehydroxyfluorene_tne2','twohydroxyfluorene_tne2','twoCyEMA_TNE2', 'allGroup')]
finalModel5 <- finalModel[,c('LBXCOT','LBXHCT','urinaryTNE2','URXNAL','twohydroxyfluorene','threehydroxyfluorene', 'twoCyEMA','NNAL_TNE2','threehydroxyfluorene_tne2','twohydroxyfluorene_tne2','twoCyEMA_TNE2', 'allGroup')]
finalModel6 <- finalModel[,c('NNAL_TNE2','threehydroxyfluorene_tne2','twohydroxyfluorene_tne2','twoCyEMA_TNE2', 'allGroup')]

set.seed(40)
finalModeltreeResult <- randomForest(allGroup ~ .,
                                 data = finalModel,
                                 ntree = 500,
                                 mtry = 16, # 自变量最大数
                                 importance = T)
# 下面每个预测
set.seed(40)
finalModeltreeResult1 <- randomForest(allGroup ~ .,
                                     data = finalModel1,
                                     ntree = 500,
                                     mtry = 5, # 自变量最大数
                                     importance = T)
set.seed(40)
finalModeltreeResult2 <- randomForest(allGroup ~ .,
                                     data = finalModel2,
                                     ntree = 500,
                                     mtry = 4, # 自变量最大数
                                     importance = T)
set.seed(40)
finalModeltreeResult3 <- randomForest(allGroup ~ .,
                                     data = finalModel3,
                                     ntree = 500,
                                     mtry = 7, # 自变量最大数
                                     importance = T)
set.seed(40)
finalModeltreeResult4 <- randomForest(allGroup ~ .,
                                     data = finalModel4,
                                     ntree = 500,
                                     mtry = 7, # 自变量最大数
                                     importance = T)
set.seed(40)
finalModeltreeResult5 <- randomForest(allGroup ~ .,
                                     data = finalModel5,
                                     ntree = 500,
                                     mtry = 11, # 自变量最大数
                                     importance = T)
set.seed(40)
finalModeltreeResult6 <- randomForest(allGroup ~ .,
                                      data = finalModel6,
                                      ntree = 500,
                                      mtry = 4, # 自变量最大数
                                      importance = T)

#可以认为直接拿着  OOB estimate of  error  当做overall 剩下那几个class error 当做误差 也就是 准确度进行计算
board <- read.csv('D:/Users/GCR/Desktop/bootBoard.csv')

# 进行拼接表格数据 2023年5月23日13:44:21
conf <- finalModeltreeResult$confusion[,-ncol(finalModeltreeResult$confusion)]
oob <-  round((sum(diag(conf))/sum(conf)),2)
oveall <-  oob * 100
numva<-finalModeltreeResult$mtry
as <-as.data.frame(finalModeltreeResult$confusion)
board[1,2] <-numva
board[1,3] <-oveall
board[1,4] <-(1-as['NEG','class.error'])*100
board[1,5] <-(1-as['TEG','class.error'])*100
board[1,6] <-(1-as['MEG','class.error'])*100
# 完成一个
conf <- finalModeltreeResult1$confusion[,-ncol(finalModeltreeResult1$confusion)]
oob <-  round((sum(diag(conf))/sum(conf)),2)
oveall <-  oob * 100
numva<-finalModeltreeResult1$mtry
as <-as.data.frame(finalModeltreeResult1$confusion)
board[2,2] <-numva
board[2,3] <-oveall
board[2,4] <-(1-as['NEG','class.error'])*100
board[2,5] <-(1-as['TEG','class.error'])*100
board[2,6] <-(1-as['MEG','class.error'])*100
# 完成2个
conf <- finalModeltreeResult2$confusion[,-ncol(finalModeltreeResult2$confusion)]
oob <-  round((sum(diag(conf))/sum(conf)),2)
oveall <-  oob * 100
numva<-finalModeltreeResult2$mtry
as <-as.data.frame(finalModeltreeResult2$confusion)
board[3,2] <-numva
board[3,3] <-oveall
board[3,4] <-(1-as['NEG','class.error'])*100
board[3,5] <-(1-as['TEG','class.error'])*100
board[3,6] <-(1-as['MEG','class.error'])*100
# 完成3个
conf <- finalModeltreeResult3$confusion[,-ncol(finalModeltreeResult3$confusion)]
oob <-  round((sum(diag(conf))/sum(conf)),2)
oveall <-  oob * 100
numva<-finalModeltreeResult3$mtry
as <-as.data.frame(finalModeltreeResult3$confusion)
board[4,2] <-numva
board[4,3] <-oveall
board[4,4] <-(1-as['NEG','class.error'])*100
board[4,5] <-(1-as['TEG','class.error'])*100
board[4,6] <-(1-as['MEG','class.error'])*100
# 完成4个
conf <- finalModeltreeResult4$confusion[,-ncol(finalModeltreeResult4$confusion)]
oob <-  round((sum(diag(conf))/sum(conf)),2)
oveall <-  oob * 100
numva<-finalModeltreeResult4$mtry
as <-as.data.frame(finalModeltreeResult4$confusion)
board[5,2] <-numva
board[5,3] <-oveall
board[5,4] <-(1-as['NEG','class.error'])*100
board[5,5] <-(1-as['TEG','class.error'])*100
board[5,6] <-(1-as['MEG','class.error'])*100
# 完成5个
conf <- finalModeltreeResult5$confusion[,-ncol(finalModeltreeResult5$confusion)]
oob <-  round((sum(diag(conf))/sum(conf)),2)
oveall <-  oob * 100
numva<-finalModeltreeResult5$mtry
as <-as.data.frame(finalModeltreeResult5$confusion)
board[6,2] <-numva
board[6,3] <-oveall
board[6,4] <-(1-as['NEG','class.error'])*100
board[6,5] <-(1-as['TEG','class.error'])*100
board[6,6] <-(1-as['MEG','class.error'])*100
# 完成6个
conf <- finalModeltreeResult6$confusion[,-ncol(finalModeltreeResult6$confusion)]
oob <-  round((sum(diag(conf))/sum(conf)),2)
oveall <-  oob * 100
numva<-finalModeltreeResult6$mtry
as <-as.data.frame(finalModeltreeResult6$confusion)
board[7,2] <-numva
board[7,3] <-oveall
board[7,4] <-(1-as['NEG','class.error'])*100
board[7,5] <-(1-as['TEG','class.error'])*100
board[7,6] <-(1-as['MEG','class.error'])*100
# write.csv(board,'/home/jar/zxjar/paper_helper/userdata/0/85/boardOut.csv')


# 完成7个
# tab 4表格完成
importance(finalModeltreeResult)

# jueq <-importance(finalModeltreeResult)
# jueq[1,][1]
# table(jueq)
# 
# colnames(jueq)
# matrix1 <-jueq[1,1]
# jueqdataf<-as.data.frame(jueq)
# colnames(jueqdataf)
# 
# sort(jueqdataf,)jueqdataf
# 
# 
# aimportant <- varImpPlot(finalModeltreeResult,main = "variable Importance score" ,type= 2,sort =T )
# aa <- as.data.frame(aimportant)
# 
# row.names(jueqdataf)
# finalModeltreeResult
# 解决网站 为了保存重要变量生成图片地址 https://data-flair.training/blogs/random-forest-in-r/
important <- importance(finalModeltreeResult, type=2 )
Important_Features <- data.frame(Feature = row.names(important), Importance = important[, 1])
#Cairo::CairoTIFF(file="test.tiff", width=8, height=8,units="in",dpi=150)
plot_ <-ggplot(Important_Features,
                aes(x= reorder(Feature,Importance) , y = Importance) ) +
  geom_bar(stat = "identity",
           fill = "#800080") +
  coord_flip() +
  theme_light(base_size = 20) +
  xlab("Questionnaire and Biomarker Variables") +
  ylab("variable Importance score")+
  ggtitle("variable Importance score") +
  theme(plot.title = element_text(size=18))

 ggsave("important_features.tiff", plot_)


# 最后一个表格复现  只是区分MEG 与 TEG

finalModel5 <- agerange317dropna[,c( 'smokersize','LBXCOT','LBXHCT','URXNAL','allGroup')]
colnames(finalModel5)

colnames(finalModel5)[1:4] <- c('smokersize','SerumCotinine','serumHydroxycotinine','UrinaryNNAL')

 finalModel5$allGroupint <- ifelse(finalModel5$allGroup=='MEG',1,0)
 finalModel5 <- missRanger(finalModel5, pmm.k = 2, num.trees = 500)
finalModel5 <- finalModel5[which(finalModel5$allGroup!='NEG'),]
finalModel5 <- finalModel5[which(finalModel5$smokersize!='0 smokers'),]
# table(finalModel5$allGroup)
# finalModel5[which(finalModel5$allGroup!='NEG'),]
# dim(finalModel5)
# dim(finalModel5[which(!is.na(finalModel5$allGroup)),])
# dim(finalModel5[which(finalModel5$allGroup!='MEG'),]) + dim(finalModel5[which(finalModel5$allGroup!='TEG'),])

# aaa <- finalModel5[which(finalModel5$allGroup!='NEG'),]


# smokersizes <- recode_factor(finalModel5$smokersize, 
#                             '0 smokers' = 0,
#                             `1 smokers` = 1,
#                             `≥2 smokers` = 2
# )
# finalModel5$smokersize <-smokersizes
# finalModel5$smokersize <-as.numeric(smokersizes)

# smokersize','SerumCotinine','serumHydroxycotinine','UrinaryNNAL
m1 <- glm(
  allGroupint ~ SerumCotinine + serumHydroxycotinine+UrinaryNNAL+smokersize, 
  data = finalModel5,
  family = binomial(link = 'logit') 
    
)

reg <- tbl_regression(m1,exponentiate =T )
reg
fit.result<-summary(m1)
df1<-fit.result$coefficients
df2<-confint(m1)
df3<-cbind(df1,df2)
df4<-data.frame(df3[-1,c(1,4,5,6)])
df4$Var<-rownames(df4)
colnames(df4)<-c("OR","Pvalue","OR_1","OR_2","Var")
df5<-df4[,c(5,1,2,3,4)]
df5$OR_mean<-df5$OR
df5$OR<-paste0(round(df5$OR,2),
               "(",
               round(df5$OR_1,2),
               "~",
               round(df5$OR_2,2),
               ")")
df5$Pvalue<-round(df5$Pvalue,3)


Cairo::CairoTIFF(file="table6.tiff", width=8, height=8,units="in",dpi=150)
a<-forestplot(labeltext=as.matrix(df5[,1:3]),
           mean=df5$OR_mean,
           lower=df5$OR_1,
           upper=df5$OR_2,
           zero=0.25,
           boxsize=0.2,
           graph.pos=2)
plot(a)









